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This  unclassified  document  reports  the  work  performed  by 
Norbert  N.  Bojarski  -  Physicist,  Consultant,  16  Circle  Drive,  Moorestown, 
New  Jersey  08057,  (609  235  3001),  under  USAF  Contract  F33615-71-C-1576, 
Fast  Fourier  Transform  Applied  to  Radar  Scattering. 
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1  April  1971  to  31  May  1972.  This  report  was  submitted  by  the  author 
on  June  1972. 
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ABSTRACT 


The  capabi 11  ty  of  prediction  of  the  radar  echoes  (radar  cross  sections) 
of  aerospace  vehicles  is  necessary  for  either  modifying  ( reducing )  such  cross 
sections  for  the  purpose  of  reducing  detectability ,  and/or  categorizing  such 
cross  sections  for  the  purpose  of  identification.  Such  predictions  car*  be 
accomplished  by  measurements  (radar  range,  unechoic  chamber,  etc.) .  The 
purpose  of  the  capability  of  computing  radar  cross  sections  ( vis-a-vis 
measurements)  is  to  reduce  both  cost  and  time.  The  state-of-the-art  of  such 
computational  methods  consists  of  computer  solving  the  scattering  integral 
equations  by  matrix  inversion  methods.  The  matrix  nature  of  such  formulations 
and  solutions  restricts  the  size  of  the  targets  for  which  radar  cross  sections 
can  be  calculated  on  even  the  largest  and  fastest  existing  computers  to  no 
more  than  the  order  of  one  wave  length,  and  renders  the  possible  solutions 
computer-time  consuming  and  costly .  The  purpose  of  the  k-space  method  is  the 
capability  of  rapid  and  cost-efficient  computer  calculations  of  radar  cross 
sections  of  aerospace  vehicles,  particularly  those  of  much  larger  size  than 
one  wave  length.  The  technical  means  by  which  this  purpose  is  achieved  is 
summarized  next. 

The  initial-value  problem  is  solved  by  means  of  a  k-space  formulation 
of  the  field  equations,  thereby  replacing  the  conventional  integral  equation 
fo emulation  by  a  set  of  two  simultaneous  algebraic  equations  in  two  unknowns 
in  two  spaces  (the  constitutive  boundary  condition  being  an  algebraic  equation 
in  x-space) .  These  equations  are  solved  by  an  iterative  generalized- 
relaxation  method  with  the  aid  of  the  Past  Fourier  Transform  (FFT)  algorithm 
connecting  the  two  spaces,  requiring  trivial  initial  approximations.  Since 
algebraic  and  FFT  equations  are  used,  the  number  of  arithmetic  multiply-add 
operations  and  storage  allocations  required  for  a  numerical  solution  are 
reduced  from  the  order  of  and  respectively  (for  solving  the  matrix 
equations  resulting  from  the  conventional  integral  equations)  to  the  order 
of  NJog2N  and  N  respectively  (where  N  is  the  number  of  data  points  required 
for  the  specification  of  the  problem) .  The  convergence  rate  of  the  iterative 
process  is  optimized  by  generalizing  the  conventional  relaxation  factors  to 
a  relaxation  function  and/or  its  generalized  inverse  determined  by  the  Eigen 
values  of  the  appropriate  Green's  function.  These  Eigen  values  are  obtained 
numerically  in  Nlog2N  operations  (vis-a-vis  the  conventionally  required  N 
operations)  by  means  of  the  FFT  of  the  Green's  function  cast  into  a  circulant 
matrix  form.  The  advantage  gained  in  speed  and  storage  is  thus  of  the  order 
of  H2/log'> N  and  N  respectively.  This  method  is  thus  considerably  more  effi¬ 
cient,  and  permits  exact  numerical  solutions  for  much  larger  problems,  than 
is  possible  with  the  conventional  integral  equation  -  matrix  inversion  method. 
Arguments  are  presented  towards  the  view  that  the  field  equations  are  more 
fundamental  in  k-space.  The  physical  and  mathematical  meaning  of  both  the 
continuous  and  discrete  k-space  representations  are  discussed.  The  details 
and  .some  numerical  results  of  the  application  of  this  method  to  three-dimon-’ 
sional  electromagnetic  scattering  are  presented.  It  is  shown  that  the 
scattered  far  fields  are  yielded  directly  in  k-space. 
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INTRODUCTION 


The  main  subjects  of  this  report  are  a  Generalized  Relaxation  method 
and  a  Generalized  Inverse  method  ftr.  iteratively  solving  large  integral 
equations.  These  methods  are  sp.  ''if  .cally  applicable  to  the  k-space 
formulation  of  such  Integral  equations.  This  k-space  formulation  was 
developed  earlier  by  this  ;  uthor  under  a  previous  contract,  and  was 
specifically  addressed  to  the  problem  of  electromagnetic  scattering .  The 
purpose  of  the  generalized  relaxation  method  and  the  generalized  inverse 
method  is  the  acceleration  of  the  convergence  rate  and  the  enlargement  of  the 
domain  of  applicability  of  the  previously  developed  iterative  method  of 
solution  of  the  k-space  formulation.  Since  the  basic  k-space  formulation 
and  its  method  of  solution  is  applicable  to  all  the  initial  value  problems 
of  mathematical  physics  which  conventionally  lead  to  integral  equations, 
this  report  consists  of  a  brief  development  of  the  k-spaoe  formulation  of 
the  initial  value  problem,  with  brief  treatment  of  the  electromagnetic 
scattering  problem  as  a  special  case;  and  with  detailed  emphasis  on  the 
development  of  the  generalized  relaxation  function  method  and  the  generalized 
inverse  method  (the  reader  is  referred  to  the  above  cited  earlier  report  for 
a  detailed  and  tutorial  presentation  of  the  k-space  formulation  of  the 
electromagnetic  scattering  problem). 

The  organization  of  this  report  is  as  follows.  In  Section  1,  the 
conventional  matrix  formulation  of  the  ini< ial  value  problem  is  briefly 
reviewed.  In  Section  2,  the  k-space  formulation  of  the  initial  value  problem 
and  its  conventional  iterative  relaxation  method  of  solution  are  formally 


Bcjarski,  N.  N. ,  K-Space  Formulation  of  the  Electromagnetic  Scattering 
Problem ,  Air  Force  Avionics  Laboratory,  Wright-Patterson  Air  Force 
Base,  Technical  Report.  AFAL-TR-71-75,  March  1971,  Final  Report  to 
USAF  Contract  F33615-70-C-1345,  AD  882  040. 


?SBr3jBgg: 


?S*I 


developed.  In  Section  3,  the  special  case  of  the  k-space  formulation  of  the 
general  wave  scattering  problem  is  formally  developed.  Section  4  consists 
of  the  detailed  development  of  the  matrix  theory  perspective  of  the  k-space 
formulation,  necessary  for  the  subsequent  development  of  the  generalized 
relaxation  method  and  the  generalized  inverse  method.  Section  5  consists  of 
the  detailed  development  of  the  generalized  relaxation  method,  and  Section  6 
consists  of  the  detailed  development  of  the  generalized  inverse  method. 

In  Section  7,  the  special  case  of  the  k-space  formulation  of  the  electro¬ 
magnetic  scattering  problem  is  formally  developed,  and  a  brief  review  of 


earlier  obtained  numeiical  results  in  summarized.  In  Section  8,  a  summary 
of  the  state-of-development  of  the  method,  its  present  limitations,  and 
recommendations  for  future  needed  research,  are  presented.  The  applicability 
of  the  method  to  the  non-monochromatic  (wide  band)  electromagnetic  scattering 
problem,  electromagnetic  problems  other  than  scattering,  and  the  general 
initial  value  problems  of  mathematical  physics,  are  also  discussed. 

This  author  is  indebted  to  Dr.  Charles  H.  Krueger,  Jr.,  of  the 
Air  Force  Avionics  Laboratory,  Wright -Patterson  Air  Force  Base,  Ohio,  for 
developing  the  details,  programming  and  verifying,  the  two-dimensional 
k-space  formulation  of  the  electromagnetic  scattering  problem. 

This  author  is  also  indebted  to  Messrs.  Leslie  E.  Whitford  and 
William  Hopkins  of  the  Computing  and  Information  Systems  Division,  Computer 
Science  Center,  Wright-Patterson  Air  Force  Base,  Ohio,  for  directing  the 
programming  effort  and  performing  the  actual  programming  respectively,  of 
the  one-  and  three-dimensional  k-space  formulations  of  the  electromagnetic 
scattering  problem,  as  well  as  for  several  important  suggestions. 
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1 .  THE  CONVENTIONAL  MATRIX  FORMULATION  OF  THE  INITIAL  VALUE  PROBLEM 


Consider  the  n- dimensional  scalar,  vector,  or  tensor  field  <f>(x)  and 
the  source  density  w(X),  governed  by  the  linear  m-th  order  differential 
field  equation 


L»  *<x>  -  -  L,  w<x> 


where  the  n-dimensional  linear  differential  scalar,  vector,  or  tensor 


operators  [_  and  are  of  the  form 


m 

L* =  £  ai  (£:)'  '■ 

1=o  '  J  ' 

m 

i  =n  J 


subject  to  the  n-dimensional  scalar,  vector,  or  tensor  constitutive  equation 


w(x)  =  a(x)  $(x)  . 


The  conventional  n-dimensional  integral  representation  of  the 
generalized  initial  value  problem  associated  with  the  field  equation  (1)  is 


<j>(x)  =  fg(xjx')  w(x')  dnx‘  +  <f>!(x)  , 

J  n 


3 


subject  to  the  constitutive  equation  (4)  where  1  (X )  is  the  externally 
imposed  field  (initial  value  field,  the  source  distribution  of  which  is 
external  to  the  problem  defined  by  the  constitutive  equation  (4)  ),  D  is 
the  domain  of  non-vanishing  o(x),  and  g(X)  is  the  appropriate  Green's 
function,  vector,  or  tensor,  satisfying  the  differential  equation 


L  9<*>  -  -  Lw  6t*> 


The  constitutive  equation  (4)  can  also  be  viewed  as  the  boundary 
conditions  imposed  by  a  specific  physical  situation  on  the  differential 
field  equation  (1),  which  is  invariant  to  the  specific  physical  situation. 

The  conventional  numerical  method  of  solution  of  this  initial  value 
problem  is  by  means  of  numerical  matrix  inversion  methods  [1],  applied  to 
the  Fredholm  integral  equation  of  the  Second  Kind,  formed  by  combining 
(4)  and  (5),  i.e. , 


<p(x)  -  |  K( x | X * )  4>(x' )  dV  =  4>'  (x)  , 


where  the  integral  transform  kernel  K(xjx')  is  given  by 


K(x|x')  =  g(x|x')  a(x')  . 


(It  should  be  noted  that  in  cartesian  coordinates,  this  kemal  K(x|x') 
is  always  a  compound  kernel  of  the  fora  K(X | x ' )  =  g(x-x')  a(x');  i.e.t 
composed  of  a  difference  kernel  and  a  separable  produot  kernel) . 
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Such  matrix  inversion  methods  require  of  the  order  of  N2  computer 
storage  allocations,  and  of  the  order  of  N3  arithmetic  multiply-add 
operations  (for  matrix  inversion)  for  the  execution  of  a  numerical  solution, 
where  N  is  the  number  of  data  points  required  for  the  numerical  specifi¬ 
cation  of  the  constitutive  equation  (4)  (the  specification  of  the  non-van¬ 
ishing  portion  of  a(x)  ).  The  practical  size  limit  with  state-of-the-art 
computers  is  for  N  of  the  order  of  several  hundred. 
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2.  THE  K-SPACE  FORMULATION  OF  THE  GENERAL  INITIAL  VALUE  PROBLEM 


The  kspace  representation  and  solution  of  the  generalized  n- 
dimensional  initial  value  problem  is  presented  next. 

The  n-dimensional  Fourier  Transform  of  the  differential  field  equation 
(1)  yields  the  local  algebraic  scalar,  vector,  or  tensor  kspace  field 
equation 


Lw(k)  =  2  bj  Uk)J  (12) 

j=o 


Preceding  page  blank 
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The  k-space  representation  of  the  generalized  n-dimemional  initial 
value  problem 3  consistent  with  the  x-space  integral  representation  (5),  thus 
is  the  algebraic  scalar,  vector,  or  tensor  equation  (vis-a-vis  the  conven¬ 
tional  integral  or  differential  equation  representation) 


$(k)  =  G(k)  W( k)  +  $'<k)  , 


(13) 


subject  to  the  algebraic  x-space  constitutive  equation  (4),  i.e.3 


w(X)  =  a(X)  $(X)  , 


(14) 


where. 


G(k)  = 


Lw(kl 

L  (k) 
<p 


(15) 


which  clearly  can  be  taken  as  the  Green’s  function  (or  vector  or  tensor)  in 
k-space,  i.e.3 


G(k)  g(X) 


(16) 


The  generalized  n-dimensicmal  initial  value  problem  is  thus  reduced 
to  a  set  of  two  local  algebraic  ( scalar 3  vector3  or  tensor)  equations  in  two 
unknowns  in  two  spaces3  i.e.3  (13)  and  (14). 

The  k-space  representation  (13)  et  eeq.  also  clearly  follows  from  (S) 
and  the  n-dimensional  convolution  theorem  since 
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g(xjx')  =  g(X-X') 


(17) 


in  any  cartesian  coordinate  system. 

The  unique  existence  [2]  of  this  k-space  representation  is  restricted 
to  media  for  which 


(18) 


If  o(x)  is  in  general  non-vanishing  only  in  a  finite  n-dimensional 
x-domain,  then  the  pair  of  algebraic  equations  (13)  and  (14)  can  be  solved 
numerically  with  the  aid  of  the  n-dimensional  Fast  Fourier  Transform  (FFT) 
algorithm  [3]  as  the  connection  between  the  tuo  spaces,  by  the  following 
iterative  relaxation  method  [4];  the  recursion  relationship  for  which  is 


W  (k)  =  p(k|x)  w  (x) 
n  n 

(19.1) 

$  (k)  =  G(k)  W  (k)  + 
n  n 

$'(k) 

(19.2) 

«t>n (x)  =  p(x|k)  4>n (k) 

(19.3) 

n+1<X)  =  cx  O(X)  *n(X> 

+  (1-a)  w  (X)  , 
n 

(19.4) 

where  a  is  an  appropriately  chosen  relaxation  coefficient  (best  numerical 
results  to  date  were  obtained  for  a=?),  and  where  p( k [ X )  and  ptx|k)  designate 
the  Fast  Fourier  Transform  algorithm  operator  and  its  inverse  respectively. 

The  initial  approximation  w0(X)  can  be  taken  as  any  known  simply 
programmable  approximation  to  the  problem,  including  the  trivial  case 
w0(x)  =  0,  the  computer  programming  of  which  consumes  virtually  no  signifi- 
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cant,  much  needed  for  data,  core  storage  allocation  (which  is  not  the  case 
with  any  other  non-trivial  initial  approximation,  which,  at  best,  would 
reduce  by  two  the  total  number  of  iterations  required). 


In  order  to  avoid  the  numerical  difficulties  arising  from  aliasing , 
i.e the  fictitious  periodic  nature  of  the  FFT  (which  is  a  discrete  and 
finite  Fourier  transform,  vis-a-vis  the  continuous  and  infinite  Fourier 
transform  implied  by  (10)),  and  the  possible  singularities  in  the  Green's 
function  (or  vector  or  tensor)  in  k-space,  it  becomes  necessary  to  choose 
an  n-dimensional  hyper- rectangular  box  of  turice  the  size  (in  each  dimension) 
of  the  smallest  hyper- rectangular  box  in  which  the  non-vanishing  o(X)  is 
imbeddable  as  the  x-domain  for  the  FFT,  and  take  the  Green's  function,  vector, 
or  tensor  as  (in  conventional  FFT  notation  [5]) 


Nj  N2 

2  2  _1 


G(u,v, • • • ,o)  =  -  A  x 


ZEE 


•••  .  (20) 
1  2  n 


a=~2  Y='2 


where 


2ui/N. 

K,  *  e  J  ,  J*1 ,2, • ,n  (21) 

J 

Ax.  Ak.  =  2tt/N  .  (22) 

J  J  J 

Anx  *  Axi  Ax2  Axn  (23.1) 

N  =  Nj  N2  •••  Nn  ,  etc.,  (23.2) 


and  where  appropriate  use  must  be  made  of  the  n-fold  periodicity  properties 

of  the  FFT  [6]  for  both  (20),  as  well  as  the  desired  placement  of  o(x)  and 
♦ 

<*>  (x)  in  the  hyper-rectangular  FFT  x-doraain. 
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The  numerical  difficulty  arising  from  the  possible  singularity  of  the 
Green's  function  at  the  origin  of  the  x-space  (i.e.,  g(0)  =  »)  can  be 
alleviated  by  taking  advantage  of  the  appropriate  principal  value  integral 
representation  of  the  field  equation  (5),  (e.p.,  [7]),  i.e.j 


<j>(x)  =  i  w(X)  +  p|  g(  x  |  x ' )  w(x') 


’/ 

J  T) 


I 

d  x' 


(x) 


(24) 


It  thus  immediately  follows  that  g(0)  for  (20)  can  be  taken  as 


Since  algebraic  and  FFT  equations  are  used,  the  number  of  arithmetic 
multiply-add  operations  required  [8]  for  a  solution  is  reduced  from  the 
order  of  N3  (required  for  the  numerical  matrix  inversion  needed  for  solving 
the  matrix  equations  resulting  from  the  conventional  integral  equation 
representation  of  the  problem)  to  the  order  of  N  Zop2  N>  and  the  storage 
requirement  is  reduced  from  the  order  of  N2  (required  for  storing  the 
matrix  associated  with  the  Green's  function  needed  for  the  matrix  method  of 
solution)  to  the  order  of  N  (required  for  storing  the  k-  and  x-space  functions) 
The  advantage  gained  in  speed  and  storage  is  thus  of  the  order  of  N2  logi  N 
and  N  respectively.  This  method  is  thus  considerably  more  efficient,  and 
permits  exact  numerical  solutions  for  much  larger  problems.  To  date, 
acoustic  scattering  problems  of  the  order  of  N= 1 04  have  been  successfully 
solved,  and  problems  of  the  order  of  107  are  feasible  with  state-of-the-art 
computers. 
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3.  THE  K- SPACE  FORMULATION  OF  THE  GENERAL  WAVE  SCATTERING  PROBLEM 


For  the  general  n-dimensional  wave  scattering  problem,  the  (range  and 
phase  normalized)  scattered  far  fields  in  x-space  are  in  general  simply  and 
algebraically  related  to  the  n-dimensional  Fourier  Transform  of  the  induced 
(by  the  incident  field  4>'(X)  )  source  distribution  w'X),  i.e W(k),  which 
is  clearly  yielded  directly  by  the  iterative  solution  (19)  without  addi¬ 
tional  computations.  Since  this  is  not  the  case  with  the  conventional 
matrix  method  of  solution  of  the  integral  equation  representation  of  the 
scattering  problem,  the  k-space  method  of  solution  presented  is  particularly 
and  additionally  attractive  when  applied  to  scattering  problems. 

For  the  special  case  of  the  n-dimensional  Helmholtz  (time-reduced) 
wave  equations  for  which 


^i)(%)2+k5-  t26) 

j=o  J 

where  k0  5  ^  ,  and  c  is  the  wave  velocity  in  the  free  space  (this  deviation 
from  conventional  notation  is  for  the  purpose  of  distinction  from  k,  the 
Fourier  Transform  variable  of  x) ,  the  n-dimensional  k-space  Green's  function, 
in  the  notation  of  (9),  (15),  and  (16),  is  clearly 


G(k)  =  'J'(k)  Lw(k>  ,  (27.1) 

1 

¥<k)  =  -  ,  (27.2) 

k2  -  k2 

The  form  of  the  Green's  function  Y(k)  in  k-space  is  clearly  invariant 
to  the  dimensionality  n  of  the  space,  which  is  not  the  case  for  the  Green's 
function  iji(x)  in  x-space. 


Preceding  page  blank 
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Thus,  by  virtue  of  the  previously  stated  relationship  between  the 
(range  and  phase  normalized)  scattered  far  fields  in  x-space  and  the  source 
distribution  W(k)  in  k-space,  and  conservation  of  energy  considerations 
(mathematically  equivalent  to  Parseval's  formula)  for  passive  media  of  finite 
spatial  extent  (see  (18)  ),  it  follows  that  the  Badiation  Condition  (the 
boundary  condition  or  constitutive  equations  at  infinity)  for  the  Helmholtz 
equation  Green's  function  in  k-space  can  be  stated  as 

*<IU 

- —  <  -  ,  !k  |  =  k0  ,  (28) 

¥(k  )  s 

s 

where  is  the  propagation  wave  number  vector  of  the  scattered  far  fields. 

It  can  thus  be  argued  that  the  field  equations  of  mathematical  physics 
are  more  fundamental  in  k-space  because  of  the  simple  local  algebraic  nature 
of  these  equations  in  k-space  ( vis-a-vis  the  global  nature  of  the  integral 
or  differential  representation  of  these  field  equations  in  x-space),  and 
the  invariance  of  the  form  of  the  k-space  Green's  function  for  the  Helmholtz 
equation  with  respect  to  the  dimensionality  of  the  space;  particularly  when 
bearing  in  mind  that  the  Fourier  Transform  is  the  only  transform  known  for 
which  a  fast  algorithm  exists.  However,  the  constitutive  equations  (or 
boundary  conditions)  are  more  fundamental  in  x-space  because  of  their  local 
algebraic  nature  in  x-space. 
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4.  THE  MATRIX  THEORY  PERSPECTIVE  OF  TOE  K-SPACE  FORMULATION 


The  equivalent  matrix  theory  perspective  of  the  k-space  formulation  can 
best  be  developed  from  the  integral  equation  (7)  for  the  one -dimensional 
problem,  with  the  generalization  to  n-dimensional  integral  equations  then 
becoming  obvious;  i.e.t  by  (7)  and  (8) 


4>(x)  -  I  g(x|x')  o(x')  4>(x’)  dx*  =  <j>  (x)  . 


If  x  is  discreteized  in  M  equal  intervals  Ax  in  the  domain  (o,a),  then 


x^  =  u  Ax  ;  a  =  0,1,2,*,»,M  ; 


and  (29)  yields 


n  i 

(x  )  -  g(x  lx  ,)  a(x  ,)  <|>(x  .)  Ax  =  (^(x  ) 

a  t—J  a1  a*  a'  Y  o'  a 


Introducing  the  M-dimensional  vectors  and  41^  as 


4  =  $(x  ) 
o  a 


4>T  =  (x  )  , 

a  a 
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the  M-dimensional  (square)  diagonal  matrix  a  .  as 


0 


aB 


4X  SaB  0(xa)  , 


(34) 


1  ;  a  =  e 


0  ;  a  i  6  , 


(35) 


and  the  M-dimensional  matrix  g  ,,  as 

3aS 


JaB 


=  9<xJxft>  » 


(36) 


then  in  cartesian  tensor  notation  (with  summation  convention  implied),  (31) 
becomes 


^aB  °By 


(37) 


which,  in  conventional  matrix  notation  becomes 


♦ 


go$  -  <i> 


(38) 


(I  -  ga>d> 


(39) 


where  I  is  the  identity  matrix. 
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(The  conventional  matrix  inversion  solution  of  the  problem  is 
accomplished  by  forming  a  matrix  s=gcr  directly,  without  the  separate  intro¬ 
duction  of  the  diagonal  matrix  a,  and  numerically  inverting  the  matrix  (I-s), 

_1  | 

yielding  the  solution  4>=( I— s)  <j>  ). 

Introducing  the  M-dimensional  generatrix  vector  yo  (generatrix  of 
the  matrix  g) 


yn  =  g(x  ) 
a  a. 


(40) 


and  recalling  that  the  Green's  function  g(x|x’)  is  always  a  difference 
function  (see  text  subsequent  to  8)  of  the  form 


g(x|x* )  =  g( | x— x * | )  , 


(41) 


yields  with  the  aid  of  (361 


=  Y 


a-8 


=  Y 


B-o 


=  Y 


|a-6 


(42.1) 

(42.2) 

(42. S) 


It  thus  follows  that  the  (MxM)  Green’s  matrix  g  is  a  symmetric  matrix 
derivable  from  the  M-dimensional  Green *s  generatrix  vector  y  by  (42) . 

Introducing  a  (2MH)  dimensional  space,  an  N-dimensional  space 

of  dimension  (2W-1),  i.e. 

17 

a  i  •-  = — *---••  —  -  1 —  —  _ _ _ _ _ 


N  =  2M  +  1 


pj-f; 


►^Wpw*\f!F3=3?s!t  ty ^ 


(43) 


for  the  vectors  $,  <J> ^ ,  and  y  in  which  the  first  M  components  are  defined  as 
per  (30)  through  (42) ,  and  the  remaining  components  are  defined  by  the 
periodic  relationships 


and  in  which  the  first  MxM  components 
defined  as  per  (34)  and  (35),  and  the 
are  defined  as  zero,  i.e.t 


^  ;  M  <  a  <  N-1  ,  (44) 

4*N-a  ;  M  <  a  <  N-1  ,  (45) 

Ym  ;  M  <  a  <  N-1  ,  (46) 

N-a 

of  the  diagonal  of  the  matrix  a  are 
remaining  components  of  the  diagonal 


Ax  6  o(x  )  ;  0  <  a, 3  <  M 
a6  ct 

0  ;  M  <  a,g  <  N-1  , 


(47) 


clearly  leaves  the  now  N-dimensional  matrix  equation 


$  -  go<j>  =  <J> 


(48) 


consistent  with  (31),  (38),  and  (39). 

However,  the  NxN  Green’s  matrix  g  resulting  from  the  N-dimensional 
now  periodic  generatrix  vector  y,  still  given  by  (42),  is  now  not  only 
syrmetr.f  ,  but  also  has  the  additional  property  of  being  a  circulant  matrix  [9). 
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Some  well  known  special  properties  of  circulant  matricies  [10]  will  now 
be  recast  into  a  somewhat  more  convenient  notation  and  general  fora;  for  this 
purpose,  however,  it  becomes  convenient  to  introduce  the  NxN  Fourier  matrix  F, 
defined  as 


F 

yet 


which  clearly  is  a  symmetric  unitary  matrix,  i.e. 


(49) 


F  =  F 
pa  ap 


(50) 


F  F  =  I 


(51) 


It  should  be  noted  that  this  Fourier  matrix  has  only  N  distinct 
elements  (thus  requiring  only  N  storage  allocations) .  The  matrix  multiplication 
of  this  matrix  by  a  vector  requires  only  arithmetic  multiply-add 

operations  if  the  FFT  algorithm  is  used. 

The  special  properties  of  a  circulant  matrix  can  now  be  stated  as 
follows : 

] .  The  Fourier  matrix  is  the  -unitary  transformation  matrix  which 
diagonalizes  a  circulant  matrix,  i.e . 


F  g  F+  =  g» 


(52) 


where  g*  is  the  diagonal  Eigen  matrix  of  g.  The  rationale  for  choosing 
the  factor  1/ V?T  in  (49)  is  now  evident;  namely,  the  conventional 
definition  without  this  factor  would  have  lead  to  merely  a 


similarity  transformation  matrix,  which  is  clearly  far  less  general. 


2.  The  Eigen  values  X  of  a  circulant  matrix  g  are  related  to  the 

generatrix  vector  y  of  the  incident  matrix  and  the  Fourier  matrix  by 


X 

y 


F  y 
ua  a 


(53) 


I.e.,  The  vector  formed  by  the  set  of  Eigen  values  of  a  circulant 
matrix  is  the  Fourier  transform  of  the  generatrix  vector  of  this 
circulant  matrix.  The  practical  significance  of  this  equation  (53) 
is  that  whereas  it  takes  N3  arithmetic  multiply-add  operations  to 
compute  the  complete  set  of  Eigen  values  of  a  general  matrix,  it 
takes  only  arithmetic  multiply-add  operations  to  compute 

the  complete  set  of  Eigen  values  of  a  circulant  matrix,  since  (53) 
can  be  computed  with  the  aid  of  the  FFT  algorithm  (matrix  multipli¬ 
cation  by  the  Fourier  matrix  can  always  be  accomplished  with  the 
FFT  algorithm).  Furthermore,  whereas  N2  storage  allocations  are 
needed  for  a  general  matrix,  only  N  storage  allocations  are  needed 
for  the  generatrix  vector  completely  and  uniquely  defining  a 
circulant  matrix. 

( cx) 

3.  The  Eigen  vectors  e^of  a  circulant  matrix  are  totally  independent 
of  the  (N  independent)  elements  of  the  circulant  matrix  (and  the 
elements  of  its  generatrix  vector),  and  are  proportional  to  the 
vectors  formed  by  the  rents  or  columns  of  the  Fourier  Matrix ;  i.e. 


F 


yu 


(54) 


The  practical  significance  of  this  equation  (54)  is  that  the  Eigen 
vectors  of  a  circulant  matrix  can  be  generated  directly  by  inspection 
cf  the  Fourier  matrix  (49);  i.e.,  the  N  N-th  roots  of  unity. 
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4.  The  inverse  of  a  circulant  matrix  g  is  also  a  circulant  matrix, 
and  is  given  by 


g"1  •  F+  {{}  F  , 


(55) 


where  the  quantity  ~  is  a  diagonal  matrix,  the  diagonal  elements 
of  which  consist  of  the  reciprocals  of  the  Eigen  values  \  of  the 
circulant  matrix  g. 

An  alternative  formulation  of  (55) ,  in  terms  of  the  generatrix  y 
of  the  circulant  matrix  g,  as  per  (53),  in  a  consistent  notation,  is 


-1 


(56) 


i.e.,  the  quantity  {}  in  (56)  is  a  diagonal  matrix,  the  diagonal 
elements  of  which  consist  of  the  reciprocals  of  the  elements  of 
the  vector  (v$T  Fy). 

The  practical  significance  of  (55)  and  (56)  is  similar  to  those 
stated  subsequent  to  (53). 

The  proof  of  (55)  and  (56)  follows  directly  from  (49)  et  seq. 

All  the  preceding  properties  of  a  circulant  matrix  can  be  extended  to 
n-dimensional  spaces  by  the  introduction  of  the  n-dimensional  N*N  Fourier 
matrix  (see  49)  as 


n 

n 


1 


Vj 


2tt( 

V  .a . 

Nj  J  J 


(57) 
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where  the  matrix  multiplication  of  this  n-dimensional  NxN  Fourier  matrix  by 
an  n-dimensional  vector  (of  N  elements)  is  accomplished  numerically  with  the 
aid  of  the  n-dimensional  FFT  algorithm.  (For  the  limiting  case  of  N-*»;  i.e., 
a  Hilbert  space,  the  preceding  properties  simply  revert  back  to  the  Fourier 
integral  transform  properties,  where  it  is  now  clear  that  the  Eigen  value 
spectrwr, ;  A(k)  of  a  difference  function  g(x-x')  is  the  integral  ourier  trans¬ 
form  of  the  generatrix  function  g(x)  ). 

The  iterative  solution  (19)  can  now  be  developed  in  matrix  theory 
perspective  form  (48)  et  seq .  By  (48)  and  (52) 

4  -  F+g’Fo(j)  =  4> 1  ,  (58) 


where  now  both  o  and  g’  are  diagonal,  where  the  matrix  F  has  only  N  distinct 
elements  requiring  storage  allocations,  and  where  all  the  implied  matrix 
multiplications  can  be  executed  with  the  aid  of  the  FFT  algorithm. 

Choosing  any  vector  <j>n  for  <J>  clearly  yields  the  self  consistent  vector 
4n  by  (48)  as 


<frn  =  g«4n  +  ♦ 


(59) 


If  4 

n 

taken  as  the 


is  taken  as  the  n-th  approximation  of  4,  then  clearly  can  be 
(n+1)th  approximation  of  4;  i.e. 


*n+1 


=  go$n  +  4 


(60) 


This  recursion  relationship  (60)  for  the  iterative  solution  of  (48)  is 
equivalent  to  the  Neumann  Series  expansion  of  (48) ,  the  sufficient  condition 
for  the  convergence  of  which  is  that  the  norm  of  the  matrix  operator  go  be 
less  than  unity;  i.e. 
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< 


1 


(61) 


1 9°  I 

By  (S3),  and  the  diagonality  of  a,  this  norm  can  clearly  be  computed 
in  NZo^2n  operations. 

For  those  cases  for  which  (61)  is  not  satisfied,  the  iterative 
relaxation  method  can  be  developed  as  follows;  let  a  be  a  real  scalar 
relaxation  factor  such  that 


♦n+t  =  a*n  +  (1~a)*n 


(62) 


which,  with  the  aid  of  (59)  yields 


'n+1 


=  [  I 


a(I-go)  ]4>n  +  a<J> 


(63) 


An  alternative  more  rigorous  derivation  of  (63)  is  as  follows;  (48) 
can  be  written  as 


(I 


-go )  =  d> 


(64) 


thus 


a(I-go)$ 


(65) 


$  -  $  +  a(I-go)$ 


(66) 
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(67) 


<j>  -  [  I  -  a(I-ga))4>  =  o^ 


The  Neumann  series  expansion  of  (67)  is 


=  ^  *  [  I  -  a(I-go)]J'  (otiji1 )  ; 


And,  consistent  with  (60) ,  the  recursion  relationship  for  the  iterative 
solution  of  (67)  is 


4>n+1  =11-  ct(I-go)]<)>n  +  a$ 


which  is  consistent  with  the  iterative  k-space  solution  (19).  The  sufficient 
condition  for  the  convergence  cf  (69)  is 


|  I  -  a(I-gc)  j  <  1  . 


The  conventional  relaxation  method  thus  reduces  to  finding  the 
relaxation  factor  a  which  minimizes  |I-u(I-go)|,  and  thus  maximizes  the 
convergence  rate  of  (69) . 


The  essence  of  the  k-space  formulation,  from  a  matrix  theory- 
perspective,  is  in  (52),  (58),  and  (63),  in  which  only  the  diagonal  matrioies 
3'  and  o  need  be  stored,  and  all  matrix  multiplications  (executed  by  the  FFT 
algorithm)  aTe  via  the  Fourier  matrix  F.  This  becomes  most  evident  if  (63) 
and  (52)  are  combined  into 
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I 


=  [  I  -  cxd-F^g'Fo)]^  + 


which  is  consistent  with  the  iterative  k-space  solution  (19). 


The  mathematical  meaning  [vis-a-vis  the  physical  meaning  discussed 
earlier)  of  the  k-space  formulation  is  thus  in  having  found  the  trans¬ 
formation  matrix  (the  Fourier  matrix)  which  diagonalizes  the  discrete  field 
equations  (5)  in  a  new  coordinate  system  (the  k-space),  and  possessing  the 
means  of  executing  this  transformation  economically  in  speed  and  storage 
( i.e.j  with  the  FFT  algorithm). 
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5.  THE  GENERALIZED  RELAXATION  METHOD 


The  difficulty  with  the  conventional  relaxation  method  is  the 
difficulty  of  finding  a  relaxation  factor  which  optimizes  convergence  rate, 
as  well  as  the  fasic  question  of  the  existence  of  such  a  factor  which 
assures  convergence  at  all. 

These  difficulties  are  alleviated  if  the  relaxation  factor  is  not 
restricted  to  being  a  veal  scalar,  but  is  generalized  to  being  a  complex 
matrix. 


The  trivial  and  naive  such  relaxation  matrix  a  clearly  is 


a  =  (I-ga) 


(72) 


for  which  the  iterative  recursion  relationship  (69)  and  the  convergence 
condition  (70)  respectively  become 


pn+l 


(73) 


|  I  -  a<I-go)  |  =  0  <  1 


(74) 


namely;  the  iterative  process  converges  in  one  iteration  to  the  exact 
solution.  This  trivial  and  naive  choice  of  a  as  the  solution  to  the  problem 
was  presented  merely  for  the  purpose  of  making  the  following  argument: 

Preceding  page  blank 
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The  best  knotJn  approximation  to  the  solution  of  the  problem  should  be 
ohoscn  as  a  relaxation  matrix ,  and  not  as  the  initial  approximation  t  in  the 
iterative  process  (since  convergence  is  determined  solely  by  the  relaxation 
matrix  and  is  independent  of  the  initial  approximation  chosen ). 

However,  prior  to  searching  for  a  practically  suitable  choice  of  the 
matrix  a,  certain  practical  limitations  must  be  imposed  on  the  properties  of 
this  matrix.  A  completely  arbitrary  matrix  a  will  negate  all  the  storage  and 
speed  advantages  of  the  k-space  method,  since  such  an  arbitrary  matrix  will 
require  of  the  order  of  N2  storage  allocations  and  arithmetic  multiply-add 
operations  for  the  execution  of  the  implied  matrix  multiplications  in  (69). 

A  natural  such  choice  of  special  properties  for  the  matrix  a,  which 
will  preserve  the  speed  and  storage  advantages  of  the  k-space  method,  con¬ 
sistent  with  the  matrix  theory  perspective  of  the  k-space  formulation,  clearly 
are  that  the  matrix  a  be  diagonal ,  oiraulant ,  or  compound  (the  product  of 
a  diagonal  and  a  circulant  matrix) . 

(The  conventional  scalar  relaxation  method  can  now  be  viewed  as  the 
special  case  of  the  generalized  relaxation  matrix  being  fully  degenerate) . 

The  general  required  properties  of  the  relaxation  matrix  will  be 
examined  next. 

The  norm  implied  in  the  convergence  condition  (74)  must  clearly  be 
taken  as  the  Euclidean  norm  or  the  spectral  radius t  i.e the  largest 
magnitude  (absolute  value)  of  the  Eigen  value  of  the  matrix  operator  A  in 
(74)  where 


A  =  I  -  o(I-go) 


(75) 


It  now  becomes  convenient  to  introduce  the  matrix  B,  defined  as 


I 

I 


I 

I 

1 

1 

I 
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I 

I 

I 

I 

I 

I 

1 

I 

I 

l 

I 

l 

I 

I 

I 

I 

I 

I 

t 


B  =  go  ; 


(76) 


thus 


A  =  I  -  o(I-B) 


The  Eigen  values  of  A  are  given  by  the  determinental  equation 


det  (  A  -  XA  I  )  =  0  . 


The  convergence  condition  (74)  for  optimal  convergence  rate  thus  is 


max  |*  |  <  1 


If  the  set  of  all  Eigen  values  {X^}  is  chosen  as  exactly  zero,  then 
clearly  (74)  is  satisfied  in  such  a  fashion  that  the  resulting  a  yields  an 
exact  solution  at  the  first  iteration  (as  was  the  case  with  the  previously 
presented  naive  example  of  a=(I-ga)  ^  );  i.e.t  a  numerically  closed  form 
solution  of  the  form  (see  71) 


(77) 


(78) 


(79) 


$  =  a<f> 


(80) 


For  this  choice,  the  determinental  equation  (78)  thus  becomes 
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det  A  =  0  , 


which,  with  the  aid  of  (77)  becomes 


det  [  I  -  a(I-B)  ]  =  0  , 


(81) 


(82) 


which  can  be  rewritten  as 


det  (I-o+aB)  =  0 


det  [a(a  ^-I+B)]  =  0 


det  a  •  det  [  B  -  (I-a  ^)  ]  =  0  , 


(83) 


which  is  satisfied  by 


det  [  B  -  (I-a-1)  ]  =  0 


(84) 


The  determinental  equation  which  determines  the  Eigen  values  of  the 
matrix  B  is 


det  (  B  -  Xg  I  )  =  0  , 


(85) 
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It  thus  follows  from  (84)  and  (85)  that  if  the  matrix  (I-a  )  is 
chosen  as  the  diagonal  matrix  {Ag},  i>e.t  consistent  with  the  previously 
introduced  notation  {}, 


I  -  a  ^  =  (Xg) 


then  (81)  is  satisfied  identically.  It  thus  follows  from  (86)  that  a  is  the 
diagonal  matrix 


Examination  of  (76)  thus  reveals  that  knowledge  (or  a  rapid  means  of 
computing)  of  the  Eigen  values  of  the  matrix  go  would  yield  a  numerically  closed 
form  solution  of  the  problem;  i.e.t  (80).  Such  knowledge  or  means  are  clearly 
not  available.  However,  since  o  is  a  diagonal  matrix,  an  approximation  for 
the  Eigen  values  of  the  matrix  go  are  available;  namely,  the  product  of  the 
Eigen  values  of  the  matricies  g  and  a  (which,  in  general,  for  non-diagonal  o, 
are  not  equal  to  the  Eigen  values  of  the  product  of  the  matricies),  i.e.t 


XD  -  A  A 
B  go 


where  the  Eigen  values  A^  are  the  diagonal  elements  of  the  diagonal  matrix  o, 
and  the  Eigen  values  A^  are  given  by  (53),  and  are  computable  easily  and 
rapidly.  To  the  extent  to  which  (88)  is  a  good  approximation,  a  good  choice 


for  the  diagonal  matrix  a  is  thus 


1  1  -  A  A 
\  g  a 
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It  can  be  shown  that  for  a  chosen  as  per  (89) ,  the  norm  of  A  is 
always  less  than  unity. 

In  the  k-space  notation  of  (19),  with  the  aid  of  (53),  (89)  yields 


a(x)  = 


1 

1  -  o  (x)  G(x) 


(90) 


where  G(x)  is  G(k)  colocated  in  x-space;  i.e.,  G(k)  evaluated  at  k=x. 

For  purposes  of  computer  programming  (storage  allocations)  the 

s 

solution  (19)- (90)  can  be  simplified  by  introducing  <}>  defined  by 


(91) 


W  (k)  =  p(k|x)  w  (x)  (93.1) 

n  1  1  n 

$s(k)  =  G(k)  W  (k)  (93.2) 

n  n 

4>s ( x )  =  p(x|k)  4>S(k)  (93.3) 

n  1  1  n 

w  Xl(x)  =  8(x)  [  <{>J(x)  +  4> 1  (x)  -  G(x)  w  (x)  ]  .  (93.4) 

n+i  n  n 
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An  artificial  test  case  of  N=8  was  programmed,  with  the  following 
numerical  results.  For  cases  for  which  the  conventional  relaxation  method 
succeeded,  convergence  was  accelerated  from  about  30  iterations  by  the 
conventional  method  to  about  6  iterations  by  the  generalized  method.  For 
cases  for  which  the  conventional  relaxation  method  failed,  convergence  was 
reached  in  about  10  iterations  by  the  generalized  method  for  those  cases  for 
which  the  magnitude  of  the  largest  element  of  a  was  less  than  about  103,  and 
the  generalized  method  failed  for  those  cases  for  which  the  magnitude  of  the 
largest  element  of  a  approached  about  10 3 . 

It  thus  seems  that  the  generalized  relaxation  function  method  cannot 
overcome  the  difficulties  presented  by  a  norm  of  A  of  about  103  and  larger 
( vis-a-vis  a  norm  A  £  1  for  the  conventional  method) .  Although  higher 
precision  in  the  computations  could  undoubtedly  overcome  this  difficulty, 
such  higher  computational  precision  is  totally  undesirable  since  it  defeats 
the  very  purpose  of  the  basic  method;  i.e.,  economy  of  speed  and  storage. 

It  is  to  the  alleviation  of  this  difficulty  that  the  next  section  is  addressed. 
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6.  THE  GENERALIZED  INVERSE  METHOD 


The  generalized  relaxation  function  method  can  now  be  put  in  the 
following  perspective.  The  Neumann  series  solution  (or  the  method  of 
successive  approximations)  of  (48),  i.e.t 

4>  -  go<|>  =  <f>*  (94) 

is  applicable  if 

norm  (ga)  <  1  .  (95) 

The  conventional  matrix  method  solution  of  (48)  is  applicable  if 

norm  (go)  £  1  ,  (96) 

and  possesses  the  additional  property  of  accelerating  (over  the  Neumann  series 
solution)  and  uniforming  convergence  if  (96)  is  satisfied. 

The  generalized  relaxation  method  of  solution  of  (48)  is  applicable  if 

norm  (go)  £  IQ3  ,  (97) 

and  further  accelerates  uniform  convergence  if  (97)  is  satisfied. 
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This  section  is  addressed  to  the  problem  of  solving  (84)  subject  to 


the  condition 


norm  (go)  >  1 


For  that  purpose  let  a  be  chosen  as  the  compound  matrix 


a  =  (-go) 


i.e.j  the  negative  of  the  generalized  inverse  [11]  of  the  product  of  the 
circulant  matrix  g  and  the  diagonal  matrix  o. 

The  recursion  relationship  for  the  iterative  solution  (69)  thus  becomes 


*n+1  =  (go)*'15  (  *n  -  <frl  )  , 


(100) 


subject  to  the  norm  condition  (70),  which,  due  to  the  choice  (99)  of  a 


oe comes 


I  (~go)("n  |  <  1 


(101) 


go  |  >  l  , 


(102) 


which  is  the  desired  condition  (98) . 
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Equation  (100)  reduces  to 


pn+1 


(-1)  (-1) 
0  9 


*n  “ 


♦  ) 


(103) 


Since  the  determinant  of  g  is  in  general  non- zero,  the  inverse  of  g 
can  be  computed  economically  and  efficiently  (with  N  storage  allocations  and 
5N'iog2N  operations)  by  (55)  j  the  generalized  inverse  of  g  can  thus  be  taken 
as  the  conventional  inverse  g  ^ .  Since  a  is  a  diagonal  matrix  with  one-half 
of  its  diagonal  elements  consisting  of  zeros,  the  determinant  of  o  vanishes; 
it  is  thus  necessary  to  take  its  generalized  inverse.  A  suitable  such 
generalized  inverse  is  clearly  (in  the  previously  introduced  notation)  the 
diagonal  matrix  fj};  i.e.t  a  diagonal  matrix,  the  diagonal  elements  of  which 
consist  of  the  reciprocals  of  the  elements  of  the  diagonal  matrix  a  for  those 
values  for  which  the  diagonal  elements  of  a  do  not  vanish,  and  zero  otherwise 
(thus  yieldind  a  diagonal  generalized  inverse  with  again  one-half  of  its 
diagonal  elements  consisting  of  zeros) .  Equation  (103)  can  thus  be  written 
as 


pn+1 


-  -1  / 
=  {-}  g  ( 


-  ) 


(104) 


Close  examination  of  the  derivation  of  (48)  from  (38)  reveals  that 
(105)  could  also  have  been  derived  directly  from  (38)  by  a  process  similar 
to  the  one  that  lead  to  (48). 

Convergence  of  (104)  under  conditions  (102)  can  be  further  accelerated 
by  the  re-application  of  the  method  of  Sect.  4;  i.e.* 


♦n  *  <?>  s"' 


^n+1  c  a4>n  +  (1-a)$n 


(105.1) 

(105.2) 
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For  computer  programming  purposes,  (105)  can  clearly  be  put  into  the 
simplified  concise  form  of  (93). 

With  the  results  of  this  and  the  previous  sections,  it  is  thus  always 
possible  to  determine  the  norm  of  (go)  economically  and  efficiently,  and 
choose  the  appropriate  economical  and  efficient  method  of  solution  of  (48), 
no  matter  what  that  norm  is. 


7.  K-SPACE  FORMULATION  OF  TOE  ELECTROMAGNETIC  SCATTERING  PROBLEM 


Three-dimensional  electromagnetic  monochromatic  scattering  by  passive 
inhomogeneous  meaia,  including  perfect  conductors,  of  finite  spatial  extent 
and  arbitrary  shape,  is  considered. 


The  time-reduced  electric  and  magnetic  field  wave  equations,  valid  for 
all  linear  inhomogeneous  media,  in  terms  of  the  total  current  density  [12], 
are  respectively 


VxvxE(x)  -  k§  E(x)  =  luuo  J(x) 


(106) 


VxV>:H(X)  -  k!  H(X)  =  VxJ(x)  , 


(107) 


which,  with  the  aid  of  Maxwell's  first  and  second,  equations,  and  the  equation 
of  continuity  for  the  total  charge  and  current  density,  can  be  written  as 


=  +  ~T2  VV*J) 


(108) 


=  -VxJ  . 


(109) 


For  non-magnetic  media  and  perfectly  conducting  media,  the  appropriate 
constitutive  equations  for  the  total  volume  and  surface  current  density  0(x) 
and  K(X)  respectively  are 
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(110) 


J  =  (  af  -  luixe  >  E 


K  =  n  x  H 


(111) 


where  and  xe  are  the  free  charge  conductivity  and  electric  susceptibility 
of  a  (non -magnetic)  medium  respectively,  and  n  is  the  outward  surface  unit 
vector  of  a  perfectly  conducting  medium. 

The  latter  (111)  is  usually  regarded  as  a  boundary  condition  for 
perfect  conductors,  but  in  the  context  of  this  paper,  this  equation  must  be 
taken  as  a  constitutive  equation  in  the  truest  sense,  particularly  if 
regarded  as  a  geometrically  constraining  condition  on  the  flow  of  all  charges. 

By  the  relationship  between  surface  and  volume  current  densities, 
consistent  with  the  FFT  notation,  the  volume  current  density  for  a  perfectly 
conducting  medium,  can  be  written  as 


£  * 

dv 


(112) 


=  7-n  xH 
Av 


(113) 


=  jh  *  « 


(114) 


where  AS  is  the  finite  differential  vector  surface  area  in  the  FFT  cell  of 
volume  Av  =  A3x. 


The  conventional  magnetic  [13]  and  electric  [14]  field  principal  value 
scattering  integral  equations 
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where 


0(x) 

(117.3) 

G(k)  =  F(k|x)  Vg(X) 

(117.4) 

( 4iKr) 

7g(X)  •  |  \  4’r  / 

;  X  *  0 

;  X  =  0 

(117.5) 

r  5  |x| 

(117.6) 

and 
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^.  A1,  4  fiyj!  jy*1 


where 


E(k)  =  I’(k)  •  J(k)  +  El(k)  (118.1) 

JCX)  »  [  c,(x)  -  tux  CX)  ]  ECx)  ,  (118.2) 

T  © 


rck)  »  ptkjx)  700 


I ujuo  Cl  +  ,T2  W)g  ; 

Ko 

3tue»  *  * 


xiiO 

x  =  0 


(118.3) 

(118.4) 


Equations  (117)  and  (118)  can  now  be  numerically  solved  economically 
and  efficiently  by  the  methods  of  the  preceding  sections. 

Defining  the  range  and  phase  normalized  scattered  far-field  S(k  )  as 

b 


SCk  )  =  Mf*  Fs(x)  e',ks*x  (119) 

S  j  s 

where  F  (x)  is  any  scattered  field  satisfying  the  relationship  F  =  F  +  F  ) 
which  is  consistent  with  the  conventional  definition  [15]  of  the  radar  power 
cross  section  a  and  the  relationship 


0  = 


U20) 


readily  reveals  that  the  range  and  phase  normalized  electric  and  magnetic 

scattered  far-fields  Sg  and  $m  are  given  directly  by  the  k-space  current 

density  distributions  ^ 

I 
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Ik  7 

S_(k_)  =  — —  [  J(k>  -  kk-J(k)  ] 


e  s 


y^ir 


k=k 


i 

S  (k  )  =  —  kxj(k) 
m  5 


k=k 


(121) 


(122) 


where  z0  is  the  impedence  of  free  space.  (As  dictated  by  the  transversality 

of  the  scattered  far-fields  in  free  space.  S  ,  S  and  K  are  indeed  all 

*  e  m  s 

orthogonal  to  each  other) . 

It  can  thus  be  shown  that  the  conventionally  defined  [16]  electric 
polarization  scattering  matrix  p ^  is  given  directly  in  k-space  by 


1  ko  Zo 

5  = -  t?c*J(k) 

n5  5 


k=k 


(123) 


where,  in  the  conventional  notation  for  spherical  coordinates,  1)  are  the 
Eigen-polarizations  (spherical  coordinate  unit  base  vectors)  4>s  and  8 
associated  with  the  scattered  far- fie Id  propagation  vectc  '  kg,  and  J(k)  is 
the  k-space  current  density  induced  by  an  electric  incidt.it  plane  wave  field 
of  the  form  and  polarization 


E 1 ( x)  =  £  e,ki 


(124) 


where  £  are  the  Eigen  polarizations  (spherical  coordinate  unit  base  vectors) 

<p.  and  associated  with  the  incident  propagation  vector  k.. 

The  solution  to  the  k-space  formulation  of  the  three-dimensional 
scattering  problem  (for  the  electric  field  equations  (118)  for  non-magnetic 
media,  and  the  magnetic  field  equations  (117)  for  perfect  conductors)  has  been 
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numerically  computer  executed  for  a  limited  number  of  cases  by  the  iterative 
method  of  solution  (19),  with  final  results  within  about  one  db.  of  exact 
known  analytic  closed  form  solutions  after  about  30  iterations.  For  example, 
see  Fig.  1-5  for  comparison  of  this  technique  with  the  exact  solution  of 
Mie  [17]  for  the  perfectly  conducting  sphere  (of  radius  a) .  The  failure  of 
this  k-space  technique  in  the  near-vicinity  of  k0a  =  2.75  (see  Fig.  1)  is  due 
to  the  fact  that  k0a  =  2.75  is  the  occurrence  of  the  first  Eigen  frequency 
(internal  resonance  of  perfectly  spherical  shell).  This  difficulty  can  be 
readily  and  simply  alleviated  by  the  appropriate  incorporation  of  the  method 
of  Mitzner  [18]  into  the  k-space  method.  However,  since  the  objective  of 
this  project  was  to  prove  the  feasibility  and  merits  of  the  k-space  method, 
and  not  the  generation  of  an  operational  user-library  of  computer  programs, 
such  an  incorporation  was  taken  as  beyond  the  scope  of  this  project. 

For  the  results  of  the  application  of  the  k-space  method  to  two- 
dimensional  electromagnetic  scattering,  the  reader  is  referred  to  the  work 
of  Krueger  [19]. 
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FIGURE  1  SPHERE;  CROSS  SECTION  ve.  FREQUENCY 
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8.  SUMMARY  OF  CONCLUSIONS  AND  RECOMMENDATIONS 

The  basic  feasibility  and  merits  of  the  application  of  the  k-space 
method  to  electromagnetic  scattering  has  essentially  been  demonstrated 
sufficiently  to  warrant  the  following  conclusions  and  recommendations. 

For  small  scatterers  (N<102),  the  preference  of  the  k-space  method 
over  the  conventional  integral  equation  -  matrix  inversion  method  cannot  yet 
be  justified.  For  medium-sized  scatterers  (N~300),  the  k-space  method  is 
indeed  more  efficient  than  the  integral  equation  -  matrix  inversion  method. 

For  large  scatterers  (N>103),  the  k-space  method  is  capable  of  yielding  useful 
results  in  realistic  computer  time;  whereas  the  integral  equation  -  matrix 
inversion  method  cannot  even  be  computer  implemented  for  such  sizes. 

For  present-day,  state-of-the-art  computer  size  and  speed,  the  size 
limit  for  which  the  k-space  method  could  be  implemented  is  of  the  order  of 
N~107  {i.e.,  about  10,000  times  as  many  drta  points  than  possible  with  the 
conventional  integral  equation  ■  matrix  inversion  method) . 

For  the  generation  of  a  user-library  type  computer  program  of  the  k- 
space  method  that  would  utilize  a  maximum  of  the  basic  inherent  advantages 
and  applicabilities  or  xne  method,  the  following  effort  would  be  needed. 

1.  Implement  appropriately  and  efficiently  the  generalized  relaxation 
function  and  the  generalized  inverse  methods  into  the  k-space 
formulation  of  the  electromagnetic  scattering  program. 

2.  Develop  a  completely  general  computer  program  for  reading 
arbitrary  shapes  and  arbitrary  electromagnetic  properties  into 
the  k-space  method. 

3.  Conduct  a  thorough  error  analysis  such  that  error  bounds  will  be 
available  to  the  user  of  such  a  system. 
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4.  Develop  (or  implement  .an  existing)  Fast  Fourier  Transform  (FFT) 
program  operative  on  disc-to-core-to-disc  data,  which  would  be 
capable  of  handling  of  the  order  of  107  data  points;  i.e.t  the 
state-of-the-art  size  limit  of  discs  ( vis-a-vis  the  presently 
utilized  in-core  FFT  limited  to  about  101*  data  points). 

5.  Investigate  the  feasibility  and  desirability  of  replacing  a  disc- 
to-core-to-disc  soft-wired  (software  compiled)  FFT  program  (see 
item  4  above)  with  a  hard-wired  (hardware  compiled)  Fast  Fourier 
Analyzer  (FFA) . 

6.  Implement  a  mixed-mode  k-space  formulation  capable  of  solving 
three  simultaneous  equations  in  three  unknowns  ( i.e.3  the  electric 
and  magnetic  fields  simultaneously  with  the  current  densities), 
applicable  simultaneously  to  volume,  surface,  and  line  current 
density  representations  of  the  electromagnetic  scattering  problem; 
thus  realizing  solutions  for  large  complex  compound  scattering 
problems  ( e.g perfect  conductors  with  absorbing  materials  and 
wire  antennas). 


The  preceding  conclusions  and  recommendations  dealt  with  monochromatic 
electromagnetic  scattering;  the  subject  to  which  this  report  was  addressed. 

It  is  now  reasonable  to  conclude  that  the  k-space  method  is  also  applicable 
with  similar  advantages  to  wide  hand  electromagnetic  scattering .  This  can 
be  accomplished  by  a  four-d'"  nensional  k-space  formulation  of  the  relativistic 
four-dimensional  time-dependent  Poisson  equation;  vis-a-vis  the  presently 
implemented  k-space  formulation  of  the  Helmholtz  equation  (time-reduced  wave 
equation).  The  additional  advantage  of  such  a  formulation  is  that  the 
relativistically  correct  Doppler  shifted  spectrum  for  arbitrary  time- 
dependent  motion  (including  time-dependent  rotation,  acceleration  and 
deformation)  would  be  yielded  directly  and  efficiently. 


Furthermore,  the  k-space  method  could  similarly  be  applied  with  full 
advantage  to  electromagnetic  problems  other  than  scattering;  e.g.,  radiation 
(antennas)  problems,  propagation  problems,  etc. 

Since  the  k-space  method  does  not  require  the  constitutive  boundary 
equation  conditions  to  be  linear  algebraic,  this  method  becomes  applicable 
with  full  advantage  to  (electromagnetic)  initial  value  problems  resulting 
from  interactive  systems  that  conventionally  yield  several  coupled 
simultaneous  non-linear  integro-differential  equations  (e.g.,  magneto- ionic 
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plasma  with  acoustic  waves  and  thermodynamic  coupling);  whereas  the 
conventional  matrix  inversion  method  is  not  applicable  to  such  non-linear 
integral  equations. 

In  conclusion,  it  is  noteworthy  that  the  k-space  formulation  is  appli¬ 
cable  with  full  advantage  to  all  the  initial  value  problems  of  mathematical 
physics  that  arise  from  linear  field  equations  subject,  to  linear  or  non¬ 
linear  constitutive  boundary  condition  equations;  which  conventionally  lead 
to  (linear  or  non-linear)  integral  equations,  including  multiple  simultaneous 
such  equations  governing  interactive  systems. 
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